{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/Qiskit/qiskit-tutorials/master/images/qiskit-heading.png\" width=\"500 px\" align=\"center\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## _*Shor's Algorithm for Integer Factorization*_ \n",
    "\n",
    "This notebook is based on an official notebook by Qiskit team, available at https://github.com/qiskit/qiskit-tutorial under the [Apache License 2.0](https://github.com/Qiskit/qiskit-tutorial/blob/master/LICENSE) license. \n",
    "Initially done by Anna Phan.\n",
    "\n",
    "In this tutorial, we first introduce the problem of [integer factorization](#factorization) and describe how [Shor's algorithm](#shorsalgorithm) solves it in detail. We then [implement](#implementation) a version of it in Qiskit.\n",
    "\n",
    "\n",
    "Your **TASK** is to follow the tutorial, as you implement Shor's algorithm (studied this week) to a real, albeit simple, integer factorization task!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integer Factorization <a id='factorization'></a>\n",
    "\n",
    "Integer factorization is the decomposition of an composite integer into a product of smaller integers, for example, the integer $100$ can be factored into $10 \\times 10$. If these factors are restricted to prime numbers, the process is called prime factorization, for example, the prime factorization of $100$ is $2 \\times 2 \\times 5 \\times 5$. \n",
    "\n",
    "When the integers are very large, no efficient classical integer factorization algorithm is known. The hardest factorization problems are semiprime numbers, the product of two prime numbers. In [2009](https://link.springer.com/chapter/10.1007/978-3-642-14623-7_18), a team of researchers factored a 232 decimal digit semiprime number (768 bits), spending the computational equivalent of more than two thousand years on a single core 2.2 GHz AMD Opteron processor with 2 GB RAM:\n",
    "```\n",
    "RSA-768  = 12301866845301177551304949583849627207728535695953347921973224521517264005 \n",
    "           07263657518745202199786469389956474942774063845925192557326303453731548268 \n",
    "           50791702612214291346167042921431160222124047927473779408066535141959745985 \n",
    "           6902143413 \n",
    "           \n",
    "         = 33478071698956898786044169848212690817704794983713768568912431388982883793 \n",
    "           878002287614711652531743087737814467999489 \n",
    "         × 36746043666799590428244633799627952632279158164343087642676032283815739666 \n",
    "           511279233373417143396810270092798736308917 \n",
    "```\n",
    "The presumed difficulty of this semiprime factorization problem underlines many encryption algorithms, such as [RSA](https://www.google.com/patents/US4405829), which is used in online credit card transactions, amongst other applications.\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shor's Algorithm <a id='shorsalgorithm'></a>\n",
    "\n",
    "Shor's algorithm, named after mathematician Peter Shor, is a polynomial time quantum algorithm for integer factorization formulated in [1994](http://epubs.siam.org/doi/10.1137/S0097539795293172). It is arguably the most dramatic example of how the paradigm of quantum computing changed our perception of which computational problems should be considered tractable, motivating the study of new quantum algorithms and efforts to design and construct quantum computers. It also has expedited research into new cryptosystems not based on integer factorization. \n",
    "\n",
    "Shor's algorithm has been experimentally realised by multiple teams for specific composite integers. The composite $15$ was first factored into $3 \\times 5$ in [2001](https://www.nature.com/nature/journal/v414/n6866/full/414883a.html) using seven NMR qubits, and has since been implemented using four photon qubits in 2007 by [two](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.99.250504) [teams](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.99.250505), three solid state qubits in [2012](https://www.nature.com/nphys/journal/v8/n10/full/nphys2385.html) and five trapped ion qubits in [2016](http://science.sciencemag.org/content/351/6277/1068). The composite $21$ has also been factored into $3 \\times 7$ in [2012](http://www.nature.com/nphoton/journal/v6/n11/full/nphoton.2012.259.html) using a photon qubit and qutrit (a three level system). Note that these experimental demonstrations rely on significant optimisations of Shor's algorithm based on apriori knowledge of the expected results. In general, [$2 + \\frac{3}{2}\\log_2N$](https://link-springer-com.virtual.anu.edu.au/chapter/10.1007/3-540-49208-9_15) qubits are needed to factor the composite integer $N$, meaning at least $1,154$ qubits would be needed to factor $RSA-768$ above.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-26T17:12:22.738849Z",
     "start_time": "2018-09-26T17:12:22.719241Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div align=\"center\"><iframe width=\"560\" height=\"315\"  align=\"centre\" src=\"https://www.youtube.com/embed/hOlOY7NyMfs?start=75&end=126\" frameborder=\"0\" allowfullscreen></iframe></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import HTML\n",
    "HTML('<div align=\"center\"><iframe width=\"560\" height=\"315\"  align=\"centre\" src=\"https://www.youtube.com/embed/hOlOY7NyMfs?start=75&end=126\" frameborder=\"0\" allowfullscreen></iframe></div>')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As Peter Shor describes in the video above from [PhysicsWorld](http://physicsworld.com/cws/article/multimedia/2015/sep/30/what-is-shors-factoring-algorithm), Shor’s algorithm is composed of three parts. The first part turns the factoring problem into a period finding problem using number theory, which can be computed on a classical computer. The second part  finds the period using the quantum Fourier transform and is responsible for the quantum speedup of the algorithm. The third part uses the period found to calculate the factors.\n",
    "\n",
    "The following sections go through the algorithm in detail, for those who just want the steps, without the lengthy explanation, refer to the [blue](#stepsone) [boxes](#stepstwo) before jumping down to the [implemention](#implemention). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### From Factorization to Period Finding\n",
    "\n",
    "The number theory that underlines Shor's algorithm relates to periodic modulo sequences. Let's have a look at an example of such a sequence. Consider the sequence of the powers of two: \n",
    "$$1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, ...$$\n",
    "Now let's look at the same sequence 'modulo 15', that is, the remainder after fifteen divides each of these powers of two:\n",
    "$$1, 2, 4, 8, 1, 2, 4, 8, 1, 2, 4, ...$$\n",
    "This is a modulo sequence that repeats every four numbers, that is, a periodic modulo sequence with a period of four.\n",
    "\n",
    "Reduction of factorization of $N$ to the problem of finding the period of an integer $x$ less than $N$ and greater than $1$ depends on the following result from number theory:\n",
    "\n",
    "> The function $\\mathcal{F}(a) = x^a \\bmod N$ is a periodic function, where $x$ is an integer coprime to $N$ and $a \\ge 0$.\n",
    "\n",
    "Note that two numbers are coprime, if the only positive integer that divides both of them is 1. This is equivalent to their greatest common divisor being 1. For example, 8 and 15 are coprime, as they don't share any common factors (other than 1). However, 9 and 15 are not coprime, since they are both divisible by 3 (and 1). \n",
    "\n",
    "> Since $\\mathcal{F}(a)$ is a periodic function, it has some period $r$. Knowing that $x^0 \\bmod N = 1$, this means that $x^r \\bmod N = 1$ since the function is periodic, and thus $r$ is just the first nonzero power where $x^r = 1 (\\bmod N)$.\n",
    "\n",
    "Given this information and through the following algebraic manipulation: \n",
    "$$ x^r \\equiv 1 \\bmod N $$\n",
    "$$ x^r = (x^{r/2})^2 \\equiv 1 \\bmod N $$\n",
    "$$ (x^{r/2})^2 - 1 \\equiv 0 \\bmod N $$\n",
    "and if $r$ is an even number:\n",
    "$$ (x^{r/2} + 1)(x^{r/2} - 1) \\equiv 0 \\bmod N $$\n",
    "\n",
    "From this, the product $(x^{r/2} + 1)(x^{r/2} - 1)$ is an integer multiple of $N$, the number to be factored. Thus, so long as $(x^{r/2} + 1)$ or $(x^{r/2} - 1)$ is not a multiple of $N$, then at least one of $(x^{r/2} + 1)$ or $(x^{r/2} - 1)$ must have a nontrivial factor in common with $N$. \n",
    "\n",
    "So computing $\\text{gcd}(x^{r/2} - 1, N)$ and $\\text{gcd}(x^{r/2} + 1, N)$ will obtain a factor of $N$, where $\\text{gcd}$ is the greatest common denominator function, which can be calculated by the polynomial time [Euclidean algorithm](https://en.wikipedia.org/wiki/Euclidean_algorithm). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Classical Steps to Shor's Algorithm\n",
    "\n",
    "Let's assume for a moment that a period finding machine exists that takes as input coprime integers $x, N$ and outputs the period of $x \\bmod N$, implemented by as a brute force search below. Let's show how to use the machine to find all prime factors of $N$ using the number theory described above. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-26T17:12:26.676176Z",
     "start_time": "2018-09-26T17:12:26.670043Z"
    }
   },
   "outputs": [],
   "source": [
    "# Brute force period finding algorithm\n",
    "def find_period_classical(x, N):\n",
    "    n = 1\n",
    "    t = x\n",
    "    while t != 1:\n",
    "        t *= x\n",
    "        t %= N\n",
    "        n += 1\n",
    "    return n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For simplicity, assume that $N$ has only two distinct prime factors: $N = pq$.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\"> <a id='stepsone'></a>\n",
    "<ol>\n",
    "<li>Pick a random integer $x$ between $1$ and $N$ and compute the greatest common divisor $\\text{gcd}(x,N)$ using Euclid's algorithm.</li>\n",
    "<li>If $x$ and $N$ have some common prime factors, $\\text{gcd}(x,N)$ will equal $p$ or $q$. Otherwise $\\text{gcd}(x,N) = 1$, meaning $x$ and $N$ are coprime. </li>\n",
    "<li>Let $r$ be the period of $x \\bmod N$ computed by the period finding machine. Repeat the above steps with different random choices of $x$ until $r$ is even.</li>\n",
    "<li>Now $p$ and $q$ can be found by computing $\\text{gcd}(x^{r/2} \\pm 1, N)$ as long as $x^{r/2} \\neq \\pm 1$.</li>\n",
    "</ol>\n",
    "</div>\n",
    "\n",
    "As an example, consider $N = 15$. Let's look at all values of $1 < x < 15$ where $x$ is coprime with $15$:\n",
    "\n",
    "|  $x$  |         $x^a \\bmod 15$       | Period $r$ |$\\text{gcd}(x^{r/2}-1,15)$|$\\text{gcd}(x^{r/2}+1,15)$ | \n",
    "|:-----:|:----------------------------:|:----------:|:------------------------:|:-------------------------:|\n",
    "|   2   | 1,2,4,8,1,2,4,8,1,2,4...     |      4     |             3            |             5             |\n",
    "|   4   | 1,4,1,4,1,4,1,4,1,4,1...     |      2     |             3            |             5             |\n",
    "|   7   | 1,7,4,13,1,7,4,13,1,7,4...   |      4     |             3            |             5             |\n",
    "|   8   | 1,8,4,2,1,8,4,2,1,8,4...     |      4     |             3            |             5             |\n",
    "|   11  | 1,11,1,11,1,11,1,11,1,11,1...|      2     |             5            |             3             |\n",
    "|   13  | 1,13,4,7,1,13,4,7,1,13,4,... |      4     |             3            |             5             |\n",
    "|   14  | 1,14,1,14,1,14,1,14,1,14,1,,,|      2     |             1            |             15            |\n",
    "\n",
    "As can be seen, any value of $x$ except $14$ will return the factors of $15$, that is, $3$ and $5$. $14$ is an example of the special case where $(x^{r/2} + 1)$ or $(x^{r/2} - 1)$ is a multiple of $N$ and thus another $x$ needs to be tried. \n",
    "\n",
    "In general, it can be shown that this special case occurs infrequently, so on average only two calls to the period finding machine are sufficient to factor $N$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "For a more interesting example, first let's find larger number N, that is semiprime that is relatively small. Using the [Sieve of Eratosthenes](https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes) [Python implementation](http://archive.oreilly.com/pub/a/python/excerpt/pythonckbk_chap1/index1.html?page=last), let's generate a list of all the prime numbers less than a thousand, randomly select two, and muliply them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-26T17:12:28.800866Z",
     "start_time": "2018-09-26T17:12:28.789601Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "semiprime N = 216773\n"
     ]
    }
   ],
   "source": [
    "import random, itertools\n",
    "\n",
    "# Sieve of Eratosthenes algorithm\n",
    "def sieve( ):\n",
    "    D = {  }\n",
    "    yield 2\n",
    "    for q in itertools.islice(itertools.count(3), 0, None, 2):\n",
    "        p = D.pop(q, None)\n",
    "        if p is None:\n",
    "            D[q*q] = q\n",
    "            yield q\n",
    "        else:\n",
    "            x = p + q\n",
    "            while x in D or not (x&1):\n",
    "                x += p\n",
    "            D[x] = p\n",
    "\n",
    "# Creates a list of prime numbers up to the given argument\n",
    "def get_primes_sieve(n):\n",
    "    return list(itertools.takewhile(lambda p: p<n, sieve()))\n",
    "\n",
    "def get_semiprime(n):\n",
    "    primes = get_primes_sieve(n)\n",
    "    l = len(primes)\n",
    "    p = primes[random.randrange(l)]\n",
    "    q = primes[random.randrange(l)]\n",
    "    return p*q\n",
    "\n",
    "N = get_semiprime(1000)\n",
    "\n",
    "print(\"semiprime N =\",N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now implement the [above steps](#stepsone) of Shor's Algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-26T17:12:29.140212Z",
     "start_time": "2018-09-26T17:12:29.130741Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "semiprime N = 216773, coprime x = 32076, period r = 107814, prime factors = 907 and 239\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "def shors_algorithm_classical(N):\n",
    "    x = random.randint(0,N) # step one\n",
    "    if(math.gcd(x,N) != 1): # step two\n",
    "        return x,0,math.gcd(x,N),N/math.gcd(x,N)\n",
    "    r = find_period_classical(x,N) # step three\n",
    "    while(r % 2 != 0):\n",
    "        r = find_period_classical(x,N)\n",
    "    p = math.gcd(x**int(r/2)+1,N) # step four, ignoring the case where (x^(r/2) +/- 1) is a multiple of N\n",
    "    q = math.gcd(x**int(r/2)-1,N)\n",
    "    return x,r,p,q\n",
    "\n",
    "x,r,p,q = shors_algorithm_classical(N)\n",
    "print(\"semiprime N = \",N,\", coprime x = \",x,\", period r = \",r,\", prime factors = \",p,\" and \",q,sep=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quantum Period Finding <a id='quantumperiodfinding'></a>\n",
    "\n",
    "Let's first describe the quantum period finding algorithm, and then go through a few of the steps in detail, before going through an example. This algorithm takes two coprime integers, $x$ and $N$, and outputs $r$, the period of $\\mathcal{F}(a) = x^a\\bmod N$.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\"><a id='stepstwo'></a>\n",
    "<ol>\n",
    "<li> Choose $T = 2^t$ such that $N^2 \\leq T \\le 2N^2$. Initialise two registers of qubits, first an argument register with $t$ qubits and second a function register with $n = log_2 N$ qubits. These registers start in the initial state:\n",
    "$$\\vert\\psi_0\\rangle = \\vert 0 \\rangle \\vert 0 \\rangle$$ </li>\n",
    "<li> Apply a Hadamard gate on each of the qubits in the argument register to yield an equally weighted superposition of all integers from $0$ to $T$:\n",
    "$$\\vert\\psi_1\\rangle = \\frac{1}{\\sqrt{T}}\\sum_{a=0}^{T-1}\\vert a \\rangle \\vert 0 \\rangle$$ </li>\n",
    "<li> Implement the modular exponentiation function $x^a \\bmod N$ on the function register, giving the state:\n",
    "$$\\vert\\psi_2\\rangle = \\frac{1}{\\sqrt{T}}\\sum_{a=0}^{T-1}\\vert a \\rangle \\vert x^a \\bmod N \\rangle$$\n",
    "This $\\vert\\psi_2\\rangle$ is highly entangled and exhibits quantum parallism, i.e. the function entangled in parallel all the 0 to $T$ input values with the corresponding values of $x^a \\bmod N$, even though the function was only executed once. </li>\n",
    "<li> Perform a quantum Fourier transform on the argument register, resulting in the state:\n",
    "$$\\vert\\psi_3\\rangle = \\frac{1}{T}\\sum_{a=0}^{T-1}\\sum_{z=0}^{T-1}e^{(2\\pi i)(az/T)}\\vert z \\rangle \\vert x^a \\bmod N \\rangle$$\n",
    "where due to the interference, only the terms $\\vert z \\rangle$ with\n",
    "$$z = qT/r $$\n",
    "have significant amplitude where $q$ is a random integer ranging from $0$ to $r-1$ and $r$ is the period of $\\mathcal{F}(a) = x^a\\bmod N$. </li>\n",
    "<li> Measure the argument register to obtain classical result $z$. With reasonable probability, the continued fraction approximation of $T / z$ will be an integer multiple of the period $r$. Euclid's algorithm can then be used to find $r$.</li>\n",
    "</ol>\n",
    "</div>\n",
    "\n",
    "Note how quantum parallelism and constructive interference have been used to detect and measure periodicity of the modular exponentiation function.  The fact that interference makes it easier to measure periodicity should not come as a big surprise. After all, physicists routinely use scattering of electromagnetic waves and interference measurements to determine periodicity of physical objects such as crystal lattices. Likewise, Shor's algorithm exploits interference to measure periodicity of arithmetic objects, a computational interferometer of sorts. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Modular Exponentiation\n",
    "\n",
    "The modular exponentiation, step 3 above, that is the evaluation of $x^a \\bmod N$ for $2^t$ values of $a$ in parallel, is the most demanding part of the algorithm. This can be performed using the following identity for the binary representation of any integer: $x = x_{t-1}2^{t-1} + ... x_12^1+x_02^0$, where $x_t$ are the binary digits of $x$. From this, it follows that:\n",
    "\n",
    "\\begin{aligned}\n",
    "x^a \\bmod N & = x^{2^{(t-1)}a_{t-1}} ... x^{2a_1}x^{a_0} \\bmod N \\\\\n",
    "& = x^{2^{(t-1)}a_{t-1}} ... [x^{2a_1}[x^{2a_0} \\bmod N] \\bmod N] ... \\bmod N \\\\\n",
    "\\end{aligned}\n",
    "\n",
    "This means that 1 is first multiplied by $x^1 \\bmod N$ if and only if $a_0 = 1$, then the result is multiplied by $x^2 \\bmod N$ if and only if $a_1 = 1$ and so forth, until finally the result is multiplied by $x^{2^{(s-1)}}\\bmod N$ if and only if $a_{t-1} = 1$. \n",
    "\n",
    "Therefore, the modular exponentiation consists of $t$ serial multiplications modulo $N$, each of them controlled by the qubit $a_t$. The values $x,x^2,...,x^{2^{(t-1)}} \\bmod N$ can be found efficiently on a classical computer by repeated squaring."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Quantum Fourier Transform\n",
    "\n",
    "The Fourier transform occurs in many different versions throughout classical computing, in areas ranging from signal processing to data compression to complexity theory. The quantum Fourier transform (QFT), step 4 above, is the quantum implementation of the discrete Fourier transform over the amplitudes of a wavefunction. \n",
    "\n",
    "The classical discrete Fourier transform acts on a vector $(x_0, ..., x_{N-1})$ and maps it to the vector $(y_0, ..., y_{N-1})$ according to the formula\n",
    "$$y_k = \\frac{1}{\\sqrt{N}}\\sum_{j=0}^{N-1}x_j\\omega_N^{jk}$$\n",
    "where $\\omega_N^{jk} = e^{2\\pi i \\frac{jk}{N}}$.\n",
    "\n",
    "Similarly, the quantum Fourier transform acts on a quantum state $\\sum_{i=0}^{N-1} x_i \\vert i \\rangle$ and maps it to the quantum state $\\sum_{i=0}^{N-1} y_i \\vert i \\rangle$ according to the formula\n",
    "$$y_k = \\frac{1}{\\sqrt{N}}\\sum_{j=0}^{N-1}x_j\\omega_N^{jk}$$\n",
    "with $\\omega_N^{jk}$ defined as above. Note that only the amplitudes of the state were affected by this transformation.\n",
    "\n",
    "This can also be expressed as the map:\n",
    "$$\\vert x \\rangle \\mapsto \\frac{1}{\\sqrt{N}}\\sum_{y=0}^{N-1}\\omega_N^{xy} \\vert y \\rangle$$\n",
    "\n",
    "Or the unitary matrix:\n",
    "$$ U_{QFT} = \\frac{1}{\\sqrt{N}} \\sum_{x=0}^{N-1} \\sum_{y=0}^{N-1} \\omega_N^{xy} \\vert y \\rangle \\langle x \\vert$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, we've actually already seen the quantum Fourier transform for when $N = 2$, it is the Hadamard operator ($H$):\n",
    "$$H = \\frac{1}{\\sqrt{2}}\\begin{bmatrix} 1 & 1 \\\\ 1 & -1 \\end{bmatrix}$$\n",
    "Suppose we have the single qubit state $\\alpha \\vert 0 \\rangle + \\beta \\vert 1 \\rangle$, if we apply the $H$ operator to this state, we obtain the new state:\n",
    "$$\\frac{1}{\\sqrt{2}}(\\alpha + \\beta) \\vert 0 \\rangle + \\frac{1}{\\sqrt{2}}(\\alpha - \\beta)  \\vert 1 \\rangle \n",
    "\\equiv \\tilde{\\alpha}\\vert 0 \\rangle + \\tilde{\\beta}\\vert 1 \\rangle$$\n",
    "Notice how the Hadamard gate performs the discrete Fourier transform for $N = 2$ on the amplitudes of the state. \n",
    "\n",
    "So what does the quantum Fourier transform look like for larger N? Let's derive a circuit for $N=2^n$, $QFT_N$ acting on the state $\\vert x \\rangle = \\vert x_1...x_n \\rangle$ where $x_1$ is the most significant bit.\n",
    "\\begin{aligned}\n",
    "QFT_N\\vert x \\rangle & = \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1}\\omega_N^{xy} \\vert y \\rangle \\\\\n",
    "& = \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} e^{2 \\pi i xy / 2^n} \\vert y \\rangle \\:\\text{since}\\: \\omega_N^{xy} = e^{2\\pi i \\frac{xy}{N}} \\:\\text{and}\\: N = 2^n\\\\\n",
    "& = \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} e^{2 \\pi i \\left(\\sum_{k=1}^n y_k/2^k\\right) x} \\vert y_1 ... y_n \\rangle \\:\\text{rewriting in fractional binary notation}\\: y = y_1...y_k, y/2^n = \\sum_{k=1}^n y_k/2^k \\\\\n",
    "& = \\frac{1}{\\sqrt{N}} \\sum_{y=0}^{N-1} \\prod_{k=0}^n e^{2 \\pi i x y_k/2^k } \\vert y_1 ... y_n \\rangle \\:\\text{after expanding the exponential of a sum to a product of exponentials} \\\\\n",
    "& = \\frac{1}{\\sqrt{N}} \\bigotimes_{k=1}^n  \\left(\\vert0\\rangle + e^{2 \\pi i x /2^k } \\vert1\\rangle \\right) \\:\\text{after rearranging the sum and products, and expanding} \\\\\n",
    "& = \\frac{1}{\\sqrt{N}} \\left(\\vert0\\rangle + e^{2 \\pi i[0.x_n]} \\vert1\\rangle\\right) \\otimes...\\otimes  \\left(\\vert0\\rangle + e^{2 \\pi i[0.x_1.x_2...x_{n-1}.x_n]} \\vert1\\rangle\\right) \\:\\text{as}\\: e^{2 \\pi i x/2^k} = e^{2 \\pi i[0.x_k...x_n]} \n",
    "\\end{aligned}\n",
    "\n",
    "This is a very useful form of the QFT for $N=2^n$ as only the last qubit depends on the the\n",
    "values of all the other input qubits and each further bit depends less and less on the input qubits. Furthermore, note that $e^{2 \\pi i.0.x_n}$ is either $+1$ or $-1$, which resembles the Hadamard transform.\n",
    "\n",
    "Before we create the circuit code for general $N=2^n$, let's look at $N=8,n=3$:\n",
    "$$QFT_8\\vert x_1x_2x_3\\rangle = \\frac{1}{\\sqrt{8}} \\left(\\vert0\\rangle + e^{2 \\pi i[0.x_3]} \\vert1\\rangle\\right) \\otimes \\left(\\vert0\\rangle + e^{2 \\pi i[0.x_2.x_3]} \\vert1\\rangle\\right) \\otimes  \\left(\\vert0\\rangle + e^{2 \\pi i[0.x_1.x_2.x_3]} \\vert1\\rangle\\right) $$\n",
    "\n",
    "The steps to creating the circuit for $\\vert y_1y_2x_3\\rangle = QFT_8\\vert x_1x_2x_3\\rangle$, remembering the [controlled phase rotation gate](../tools/quantum_gates_and_linear_algebra.ipynb\n",
    ") $CU_1$, would be:\n",
    "1. Apply a Hadamard to $\\vert x_3 \\rangle$, giving the state $\\frac{1}{\\sqrt{2}}\\left(\\vert0\\rangle + e^{2 \\pi i.0.x_3} \\vert1\\rangle\\right) = \\frac{1}{\\sqrt{2}}\\left(\\vert0\\rangle + (-1)^{x_3} \\vert1\\rangle\\right)$\n",
    "2. Apply a Hadamard to $\\vert x_2 \\rangle$, then depending on $k_3$ (before the Hadamard gate) a $CU_1(\\frac{\\pi}{2})$, giving the state $\\frac{1}{\\sqrt{2}}\\left(\\vert0\\rangle + e^{2 \\pi i[0.x_2.x_3]} \\vert1\\rangle\\right)$.\n",
    "3. Apply a Hadamard to $\\vert x_1 \\rangle$, then $CU_1(\\frac{\\pi}{2})$ depending on $k_2$, and $CU_1(\\frac{\\pi}{4})$ depending on $k_3$.\n",
    "4. Measure the bits in reverse order, that is $y_3 = x_1, y_2 = x_2, y_1 = y_3$.\n",
    "\n",
    "In Qiskit, this is:\n",
    "```python\n",
    "q3 = QuantumRegister(3, 'q3')\n",
    "c3 = ClassicalRegister(3, 'c3')\n",
    "\n",
    "qft3 = QuantumCircuit(q3, c3)\n",
    "qft3.h(q[0])\n",
    "qft3.cu1(math.pi/2.0, q3[1], q3[0])\n",
    "qft3.h(q[1])\n",
    "qft3.cu1(math.pi/4.0, q3[2], q3[0])\n",
    "qft3.cu1(math.pi/2.0, q3[2], q3[1])\n",
    "qft3.h(q[2])\n",
    "```\n",
    "\n",
    "For $N=2^n$, this can be generalised, as in the `qft` function in [tools.qi](https://github.com/Q/qiskit-terra/blob/master/qiskit/tools/qi/qi.py):\n",
    "```python\n",
    "def qft(circ, q, n):\n",
    "    \"\"\"n-qubit QFT on q in circ.\"\"\"\n",
    "    for j in range(n):\n",
    "        for k in range(j):\n",
    "            circ.cu1(math.pi/float(2**(j-k)), q[j], q[k])\n",
    "        circ.h(q[j])\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example\n",
    "\n",
    "Let's factorize $N = 21$ with coprime $x=2$, following the [above steps](#stepstwo) of the quantum period finding algorithm, which should return $r = 6$. This example follows one from [this](https://arxiv.org/abs/quant-ph/0303175) tutorial. \n",
    "\n",
    "1. Choose $T = 2^t$ such that $N^2 \\leq T \\le 2N^2$. For $N = 21$, the smallest value of $t$ is 9, meaning $T = 2^t = 512$. Initialise two registers of qubits, first an argument register with $t = 9$ qubits, and second a function register with $n = log_2 N = 5$ qubits: \n",
    "$$\\vert\\psi_0\\rangle = \\vert 0 \\rangle \\vert 0 \\rangle$$\n",
    "\n",
    "2. Apply a Hadamard gate on each of the qubits in the argument register: \n",
    "$$\\vert\\psi_1\\rangle = \\frac{1}{\\sqrt{T}}\\sum_{a=0}^{T-1}\\vert a \\rangle \\vert 0 \\rangle = \\frac{1}{\\sqrt{512}}\\sum_{a=0}^{511}\\vert a \\rangle \\vert 0 \\rangle$$\n",
    "\n",
    "3. Implement the modular exponentiation function $x^a \\bmod N$ on the function register:\n",
    "\\begin{eqnarray}\n",
    "\\vert\\psi_2\\rangle \n",
    "& = & \\frac{1}{\\sqrt{T}}\\sum_{a=0}^{T-1}\\vert a \\rangle \\vert x^a \\bmod N \\rangle\n",
    " = \\frac{1}{\\sqrt{512}}\\sum_{a=0}^{511}\\vert a \\rangle \\vert 2^a \\bmod 21 \\rangle \\\\\n",
    "& = & \\frac{1}{\\sqrt{512}} \\bigg( \\;\\; \\vert 0 \\rangle \\vert 1 \\rangle + \\vert 1 \\rangle \\vert 2 \\rangle +\n",
    "\\vert 2 \\rangle \\vert 4 \\rangle + \\vert 3 \\rangle \\vert 8 \\rangle + \\;\\; \\vert 4 \\rangle \\vert 16 \\rangle + \\;\\,\n",
    "\\vert 5 \\rangle \\vert 11 \\rangle \\, + \\\\\n",
    "& & \\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\, \\vert 6 \\rangle \\vert 1 \\rangle + \\vert 7 \\rangle \\vert 2 \\rangle + \\vert 8 \\rangle \\vert 4 \\rangle + \\vert 9 \\rangle \\vert 8 \\rangle + \\vert 10 \\rangle \\vert 16 \\rangle + \\vert 11 \\rangle \\vert 11 \\rangle \\, +\\\\\n",
    "& & \\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\, \\vert 12 \\rangle \\vert 1 \\rangle + \\ldots \\bigg)\\\\\n",
    "\\end{eqnarray}\n",
    "Notice that the above expression has the following pattern: the states of the second register of each “column” are the same. Therefore we can rearrange the terms in order to collect the second register:\n",
    "\\begin{eqnarray}\n",
    "\\vert\\psi_2\\rangle \n",
    "& = & \\frac{1}{\\sqrt{512}} \\bigg[ \\big(\\,\\vert 0 \\rangle + \\;\\vert 6 \\rangle + \\vert 12 \\rangle \\ldots + \\vert 504 \\rangle + \\vert 510 \\rangle \\big) \\, \\vert 1 \\rangle \\, + \\\\\n",
    "& & \\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\; \\big(\\,\\vert 1 \\rangle + \\;\\vert 7 \\rangle + \\vert 13 \\rangle \\ldots + \\vert 505 \\rangle + \\vert 511 \\rangle \\big) \\, \\vert 2 \\rangle \\, + \\\\\n",
    "& & \\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\; \\big(\\,\\vert 2 \\rangle + \\;\\vert 8 \\rangle + \\vert 14 \\rangle \\ldots + \\vert 506 \\rangle +  \\big) \\, \\vert 4 \\rangle \\, +  \\\\\n",
    "& & \\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\; \\big(\\,\\vert 3 \\rangle + \\;\\vert 9 \\rangle + \\vert 15 \\rangle \\ldots + \\vert 507 \\rangle +  \\big) \\, \\vert 8 \\rangle \\, +  \\\\\n",
    "& & \\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\; \\big(\\,\\vert 4 \\rangle + \\vert 10 \\rangle + \\vert 16 \\rangle \\ldots + \\vert 508 \\rangle +  \\big)  \\vert 16 \\rangle \\, +  \\\\\n",
    "& & \\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\; \\big(\\,\\vert 5 \\rangle + \\vert 11 \\rangle + \\vert 17 \\rangle \\ldots + \\vert 509 \\rangle +  \\big)  \\vert 11 \\rangle \\, \\bigg]\\\\\n",
    "\\end{eqnarray}\n",
    "\n",
    "4. To simplify following equations, we'll measure the function register before performing a quantum Fourier transform on the argument register. This will yield one of the following numbers with equal probability: $\\{1,2,4,6,8,16,11\\}$. Suppose that the result of the measurement was $2$, then:\n",
    "$$\\vert\\psi_3\\rangle = \\frac{1}{\\sqrt{86}}(\\vert 1 \\rangle + \\;\\vert 7 \\rangle + \\vert 13 \\rangle \\ldots + \\vert 505 \\rangle + \\vert 511 \\rangle)\\, \\vert 2 \\rangle $$\n",
    "It does not matter what is the result of the measurement; what matters is the periodic pattern. The period of the states of the first register is the solution to the problem and the quantum Fourier transform can reveal the value of the period.\n",
    "\n",
    "5. Perform a quantum Fourier transform on the argument register:\n",
    "$$\n",
    "\\vert\\psi_4\\rangle\n",
    " = QFT(\\vert\\psi_3\\rangle)\n",
    " = QFT(\\frac{1}{\\sqrt{86}}\\sum_{a=0}^{85}\\vert 6a+1 \\rangle)\\vert 2 \\rangle\n",
    " = \\frac{1}{\\sqrt{512}}\\sum_{j=0}^{511}\\bigg(\\big[ \\frac{1}{\\sqrt{86}}\\sum_{a=0}^{85} e^{-2 \\pi i \\frac{6ja}{512}} \\big] e^{-2\\pi i\\frac{j}{512}}\\vert j \\rangle \\bigg)\\vert 2 \\rangle\n",
    "$$\n",
    "\n",
    "6. Measure the argument register. The probability of measuring a result $j$ is:\n",
    "$$ \\rm{Probability}(j) = \\frac{1}{512 \\times 86} \\bigg\\vert \\sum_{a=0}^{85}e^{-2 \\pi i \\frac{6ja}{512}} \\bigg\\vert^2$$\n",
    "This peaks at $j=0,85,171,256,341,427$. Suppose that the result of the measement yielded $j = 85$, then using continued fraction approximation of $\\frac{512}{85}$, we obtain $r=6$, as expected. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation <a id='implementation'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-26T17:12:32.454145Z",
     "start_time": "2018-09-26T17:12:30.093921Z"
    }
   },
   "outputs": [],
   "source": [
    "from qiskit import Aer\n",
    "from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister\n",
    "from qiskit import execute, register, get_backend, compile\n",
    "from qiskit.tools.visualization import plot_histogram, circuit_drawer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned [earlier](#shorsalgorithm), many of the experimental demonstrations of Shor's algorithm rely on significant optimisations based on apriori knowledge of the expected results. We will follow the formulation in [this](http://science.sciencemag.org/content/351/6277/1068) paper, which demonstrates a reasonably scalable realisation of Shor's algorithm using $N = 15$. Below is the first figure from the paper, showing various quantum circuits, with the following caption: _Diagrams of Shor’s algorithm for factoring $N = 15$, using a generic textbook approach (**A**) compared with Kitaev’s approach (**B**) for a generic base $a$. (**C**) The actual implementation for factoring $15$ to base $11$, optimized for the corresponding single-input state. Here $q_i$ corresponds to the respective qubit in the computational register. (**D**) Kitaev’s approach to Shor’s algorithm for the bases ${2, 7, 8, 13}$. Here, the optimized map of the first multiplier is identical in all four cases, and the last multiplier is implemented with full modular multipliers, as depicted in (**E**). In all cases, the single QFT qubit is used three times, which, together with the four qubits in the computation register, totals seven effective qubits. (**E**) Circuit diagrams of the modular multipliers of the form $a \\bmod N$ for bases $a = {2, 7, 8, 11, 13}$._\n",
    "\n",
    "<img src=\"../../images/shoralgorithm.png\" alt=\"Note: In order for images to show up in this jupyter notebook you need to select File => Trusted Notebook\" width=\"500 px\" align=\"center\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we cannot run this version of Shor's algorithm on an IBM Quantum Experience device at the moment as we currently lack the ability to do measurement feedforward and qubit resetting. Thus we'll just be building the ciruits to run on the simulators for now. Based on Pinakin Padalia & Amitabh Yadav's implementation, found [here](https://github.com/amitabhyadav/Shor-Algorithm-on-IBM-Quantum-Experience)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we'll construct the $a^1 \\bmod 15$ circuits for $a = 2,7,8,11,13$ as in **E**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-26T17:12:32.467183Z",
     "start_time": "2018-09-26T17:12:32.459175Z"
    }
   },
   "outputs": [],
   "source": [
    "# qc = quantum circuit, qr = quantum register, cr = classical register, a = 2, 7, 8, 11 or 13\n",
    "def circuit_amod15(qc,qr,cr,a):\n",
    "    if a == 2:\n",
    "        qc.cswap(qr[4],qr[3],qr[2])\n",
    "        qc.cswap(qr[4],qr[2],qr[1])\n",
    "        qc.cswap(qr[4],qr[1],qr[0])\n",
    "    elif a == 7:\n",
    "        qc.cswap(qr[4],qr[1],qr[0])\n",
    "        qc.cswap(qr[4],qr[2],qr[1])\n",
    "        qc.cswap(qr[4],qr[3],qr[2])\n",
    "        qc.cx(qr[4],qr[3])\n",
    "        qc.cx(qr[4],qr[2])\n",
    "        qc.cx(qr[4],qr[1])\n",
    "        qc.cx(qr[4],qr[0])\n",
    "    elif a == 8:\n",
    "        qc.cswap(qr[4],qr[1],qr[0])\n",
    "        qc.cswap(qr[4],qr[2],qr[1])\n",
    "        qc.cswap(qr[4],qr[3],qr[2])\n",
    "    elif a == 11: # this is included for completeness\n",
    "        qc.cswap(qr[4],qr[2],qr[0])\n",
    "        qc.cswap(qr[4],qr[3],qr[1])\n",
    "        qc.cx(qr[4],qr[3])\n",
    "        qc.cx(qr[4],qr[2])\n",
    "        qc.cx(qr[4],qr[1])\n",
    "        qc.cx(qr[4],qr[0])\n",
    "    elif a == 13:\n",
    "        qc.cswap(qr[4],qr[3],qr[2])\n",
    "        qc.cswap(qr[4],qr[2],qr[1])\n",
    "        qc.cswap(qr[4],qr[1],qr[0])\n",
    "        qc.cx(qr[4],qr[3])\n",
    "        qc.cx(qr[4],qr[2])\n",
    "        qc.cx(qr[4],qr[1])\n",
    "        qc.cx(qr[4],qr[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we'll build the rest of the period finding circuit as in **D**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-26T17:12:32.481843Z",
     "start_time": "2018-09-26T17:12:32.470502Z"
    }
   },
   "outputs": [],
   "source": [
    "# qc = quantum circuit, qr = quantum register, cr = classical register, a = 2, 7, 8, 11 or 13\n",
    "def circuit_aperiod15(qc,qr,cr,a):\n",
    "    if a == 11:\n",
    "        circuit_11period15(qc,qr,cr)\n",
    "        return\n",
    "    \n",
    "    # Initialize q[0] to |1> \n",
    "    qc.x(qr[0])\n",
    "\n",
    "    # Apply a**4 mod 15\n",
    "    qc.h(qr[4])\n",
    "    #   controlled identity on the remaining 4 qubits, which is equivalent to doing nothing\n",
    "    qc.h(qr[4])\n",
    "    #   measure\n",
    "    qc.measure(qr[4],cr[0])\n",
    "    #   reinitialise q[4] to |0>\n",
    "    qc.reset(qr[4])\n",
    "\n",
    "    # Apply a**2 mod 15\n",
    "    qc.h(qr[4])\n",
    "    #   controlled unitary\n",
    "    qc.cx(qr[4],qr[2])\n",
    "    qc.cx(qr[4],qr[0])\n",
    "    #   feed forward\n",
    "    if cr[0] == 1:\n",
    "        qc.u1(math.pi/2.,qr[4])\n",
    "    qc.h(qr[4])\n",
    "    #   measure\n",
    "    qc.measure(qr[4],cr[1])\n",
    "    #   reinitialise q[4] to |0>\n",
    "    qc.reset(qr[4])\n",
    "\n",
    "    # Apply a mod 15\n",
    "    qc.h(qr[4])\n",
    "    #   controlled unitary.\n",
    "    circuit_amod15(qc,qr,cr,a)\n",
    "    #   feed forward\n",
    "    if cr[1] == 1:\n",
    "        qc.u1(math.pi/2.,qr[4])\n",
    "    if cr[0] == 1:\n",
    "        qc.u1(math.pi/4.,qr[4])\n",
    "    qc.h(qr[4])\n",
    "    #   measure\n",
    "    qc.measure(qr[4],cr[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we build the optimised circuit for $11 \\bmod 15$ as in **C**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-26T17:12:32.493545Z",
     "start_time": "2018-09-26T17:12:32.484582Z"
    }
   },
   "outputs": [],
   "source": [
    "def circuit_11period15(qc,qr,cr):\n",
    "    # Initialize q[0] to |1> \n",
    "    qc.x(qr[0])\n",
    "\n",
    "    # Apply a**4 mod 15\n",
    "    qc.h(qr[4])\n",
    "    #   controlled identity on the remaining 4 qubits, which is equivalent to doing nothing\n",
    "    qc.h(qr[4])\n",
    "    #   measure\n",
    "    qc.measure(qr[4],cr[0])\n",
    "    #   reinitialise q[4] to |0>\n",
    "    qc.reset(qr[4])\n",
    "\n",
    "    # Apply a**2 mod 15\n",
    "    qc.h(qr[4])\n",
    "    #   controlled identity on the remaining 4 qubits, which is equivalent to doing nothing\n",
    "    #   feed forward\n",
    "    if cr[0] == 1:\n",
    "        qc.u1(math.pi/2.,qr[4])\n",
    "    qc.h(qr[4])\n",
    "    #   measure\n",
    "    qc.measure(qr[4],cr[1])\n",
    "    #   reinitialise q[4] to |0>\n",
    "    qc.reset(qr[4])\n",
    "\n",
    "    # Apply 11 mod 15\n",
    "    qc.h(qr[4])\n",
    "    #   controlled unitary.\n",
    "    qc.cx(qr[4],qr[3])\n",
    "    qc.cx(qr[4],qr[1])\n",
    "    #   feed forward\n",
    "    if cr[1] == 1:\n",
    "        qc.u1(math.pi/2.,qr[4])\n",
    "    if cr[0] == 1:\n",
    "        qc.u1(math.pi/4.,qr[4])\n",
    "    qc.h(qr[4])\n",
    "    #   measure\n",
    "    qc.measure(qr[4],cr[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's build and run a circuit for $a = 7$, and plot the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-26T17:12:32.830902Z",
     "start_time": "2018-09-26T17:12:32.495963Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEaCAYAAAAG87ApAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3X2YVWW9//H3BwhM8SnFiicRBRPEQEelU5mKD/gQqKnhdUrNDPPorzopRr/SOpw8mZrZlVqSerIs8CG1sTB/GmKoqQxCGpg6IcKIKSpqpEID398fa8242cwwezF71t6z+byuay7Xvte6Z3/3Hc13r3U/KSIwMzMrVY9KB2BmZt2LE4eZmWXixGFmZpk4cZiZWSZOHGZmlokTh5mZZeLEYWZmmThxmJlZJk4cZmaWSa9KB9AVdt555xgyZEilwzAz61bmz5//SkT06+i6mkwcQ4YMoaGhodJhmJl1K5KeL+U6P6oyM7NMnDjMzCwTJw4zM8vEicPMzDJx4jAzs0ycOMzMLBMnDjMzy8SJw8zMMsktcUgaL+lpSY2SprZx/ouSnpS0UNKDkkYUnPt6Wu9pSUfmFbOZmW0sl8QhqSdwNXAUMAI4pTAxpH4VEaMiYjRwKXBFWncEMAkYCYwHrkl/n5mZVUBedxwHAI0RsSQi1gIzgYmFF0TEmwUvtwEiPZ4IzIyINRHxHNCY/j4zM6uAvNaqGgAsL3jdBBxYfJGkc4CvAr2BQwvqPlJUd0DXhGlmZh3J645DbZTFRgURV0fE7sDXgG9mqStpsqQGSQ0rV67sVLBmZta+vBJHEzCo4PVAYMUmrp8JHJelbkRMj4i6iKjr16/DVYHNzGwz5ZU45gHDJO0mqTdJZ3d94QWShhW8PAZ4Nj2uByZJ6iNpN2AY8FgOMZuZWRty6eOIiGZJ5wL3AD2BGyJikaRpQENE1APnSjoM+BewCjgtrbtI0i3AYqAZOCci1uURt5mZbSy3eRwRMSsihkfE7hFxcVp2UZo0iIgvR8TIiBgdEYdExKKCuhen9faMiLvzitmsVv3+979nzz33ZI899uCSSy7Z6PwVV1zBiBEj2GeffRg3bhzPP//u/j7Lli3jiCOOYK+99mLEiBEsXboUgI9//OOMHj2a0aNH079/f4477riNfq/VhprcAdDM2rdu3TrOOecc7r33XgYOHMj+++/PhAkTGDHi3alVY8aMoaGhga233pof//jHXHDBBdx8880AnHrqqXzjG9/g8MMPZ/Xq1fTokXz/nDt3bmv9T33qU0ycOBGrTV5yxGwL89hjj7HHHnswdOhQevfuzaRJk/jNb36zwTWHHHIIW2+9NQBjx46lqakJgMWLF9Pc3Mzhhx8OQN++fVuva/GPf/yD2bNn+46jhjlxmG1hXnjhBQYNeneg4sCBA3nhhRfavf7666/nqKOOAuCZZ55hhx124IQTTmDMmDFMmTKFdes27HK84447GDduHNttt13XfACrOCcOsy1MxEbToJDami4FN910Ew0NDUyZMgWA5uZm5s6dy+WXX868efNYsmQJP/vZzzaoM2PGDE455ZSyx23Vw4nDbAszcOBAli9/dyGHpqYm+vfvv9F19913HxdffDH19fX06dOnte6YMWMYOnQovXr14rjjjuPxxx9vrfPqq6/y2GOPccwxx3T9B7GKceIw28Lsv//+PPvsszz33HOsXbuWmTNnMmHChA2uWbBgAWeddRb19fXssssuG9RdtWoVLaszzJ49e4NO9VtvvZVjjz2WrbbaKp8PYxXhxGFVozNDRHv27Nk6FLTwj6CHiG6sV69eXHXVVRx55JHstddenHzyyYwcOZKLLrqI+vpkXu6UKVNYvXo1J5100gZt2rNnTy6//HLGjRvHqFGjiAi+8IUvtP7umTNn+jHVFkBtPe/s7urq6qKhoaHSYVgG69atY/jw4RsMEZ0xY8YG32bvv/9+DjzwwNYhonPmzGkdItq3b19Wr169yfdoGSJ66qmndulnMeuuJM2PiLqOrvMdh1WFzgwRLYWHiJqVjxOHVYXODBEFeOedd6irq2Ps2LHceeedG13vIaJm5eOZ41YVNmeI6AMPPNBatmzZMvr378+SJUs49NBDGTVqFLvvvnvr+RkzZnDmmWeWP3CzLZDvOKwqdGaIKNB67dChQzn44INZsGBB6zkPETUrLycOqwqdGSK6atUq1qxZA8Arr7zCQw895CGiZl3Ij6qsKhQOEV23bh1nnHFG6xDRuro6JkyYsMEQUYDBgwdTX1/PU089xVlnnUWPHj1Yv349U6dO3SBxzJw5k6lTp1bqo5nVHA/HNTMzwMNxzcysizhxdKAzs5kB3nzzTQYMGMC5554LwFtvvcUxxxzDhz70IUaOHOlHKGbW7ThxbELLhjd33303ixcvZsaMGSxevHiDa1o2vHniiSc48cQTueCCCzY4f+GFF/KJT3xig7Lzzz+fv/71ryxYsICHHnqIu+/2poZm1n24c3wTCmczA62zmQs7Xg855JDW47Fjx3LTTTe1vp4/fz4vvfQS48ePp6XPZeutt26t07t3b/bdd99MM6DNNseQqb+r6PsvvcRDoWuJ7zg2oTOzmdevX895553HZZdd1u71r7/+OnfddRfjxo0rX9BmZl3Mdxyb0JnZzNdccw1HH330BomnUHNzM6eccgpf+tKXWu9ozMy6AyeOTcg6m/mBBx5onc38pz/9iblz53LNNdewevVq1q5dS9++fVs72CdPnsywYcP4yle+ks+HMTMrEz+q2oTOzGb+5S9/ybJly1i6dCmXX345p556amvS+OY3v8kbb7zBlVdemevnMbP8be7IzOeff5799tuP0aNHM3LkSH7yk59sVHfChAnsvffeXf4ZijlxbEJnNrxpT1NTExdffDGLFy9m3333ZfTo0Vx33XV5fBwzy1lnRmZ+8IMf5OGHH2bhwoU8+uijXHLJJaxYsaK13u23307fvn1z/Twt/KiqA0cffTRHH330BmXTpk1rPb7vvvs6/B2nn346p59+OpA8/qrF2fqbq5KjfTzSx7paZ0Zm9u7du7V8zZo1rF+/vvX16tWrueKKK5g+fTonn3xyV3+MjeR2xyFpvKSnJTVK2mjWm6SvSlos6QlJf5C0a8G5dZIWpj/1ecVsZtYZnd1nZvny5eyzzz4MGjSIr33ta619rBdeeCHnnXde68ZmecslcUjqCVwNHAWMAE6RNKLosgVAXUTsA9wGXFpw7u2IGJ3+bPpZkJlZldickZlTpkxpLRs0aBBPPPEEjY2N3Hjjjbz00kssXLiQxsZGjj/++C6LuyN5Pao6AGiMiCUAkmYCE4HWh30RcX/B9Y8An8kpNjOzLtGZkZmF+vfvz8iRI5k7dy4rV65k/vz5DBkyhObmZl5++WUOPvhg5syZ05UfZQN5PaoaACwveN2UlrXn80DhOhxbSWqQ9IgkbxptZt1CZ0ZmNjU18fbbbwPJnjMPPfQQe+65J2effTYrVqxg6dKlPPjggwwfPjzXpAH53XG0dW/WZg+xpM8AdUDhAk+DI2KFpKHAbElPRsTfiupNBiZDsk+DmVmldXafmfPOOw9JRATnn38+o0aNqvAnSuSyH4ekjwDfjogj09dfB4iI7xZddxjwI+ATEfFyO7/rZ8BvI+K29t7P+3F0Hx5VlQ+vVWWlqLb9OOYBwyTtJqk3MAnYYHSUpDHAtcCEwqQhaUdJfdLjnYGPUtA3YmZm+crlUVVENEs6F7gH6AncEBGLJE0DGiKiHrgM6Avcmo46WJaOoNoLuFbSepJEd0lEOHGYmVVIbhMAI2IWMKuo7KKC48PaqfcwUB0P9szMzDPH2+Ln7mZm7fNaVWZmlokTh5mZZeLEYWZmmThxmJlZJk4cZmaWiUdVmZmV0ZYwS993HGZmlokTh5mZZeLEYWZmmThxmJlZJk4cZmaWiROHmZll4sRhZmaZOHGYmVkmThxmZpaJE4eZmWXixGFmZpk4cZiZWSZOHGZmlokTh5mZZeLEYWZmmZScOCT1k9Q3Pe4p6XOSTpXk5GNmtgXJ8kf/t8Cw9Phi4Hzgq8D3yx2UmZlVryw7AA4HFqbHnwH+DVgNLAL+s8xxmZlZlcpyx7EO6C1pFPBGRCwDXgf6llJZ0nhJT0tqlDS1jfNflbRY0hOS/iBp14Jzp0l6Nv05LUPMZmZWZlnuOO4GbgF2AmamZSOAFzqqKKkncDVwONAEzJNUHxGLCy5bANRFxFuSzgYuBT4t6X3At4A6IID5ad1VGWI3M7MyyXLHcSbwO+B64Ltp2c7At0uoewDQGBFLImItSeKZWHhBRNwfEW+lLx8BBqbHRwL3RsRrabK4FxifIW4zMyujku84ImINMD0dRfV+4MWImFNi9QHA8oLXTcCBm7j+8yR3OO3VHVDi+5qZWZllGY67g6RfAe8AjWnZBEnfKaV6G2XRzvt8huSx1GVZ6kqaLKlBUsPKlStLCMnMzDZHlkdVPwHeAHYF1qZlfwI+XULdJmBQweuBwIriiyQdBnwDmJDe4ZRcNyKmR0RdRNT169evhJDMzGxzZEkc44AvRcSLpN/4I2IlsEsJdecBwyTtJqk3MAmoL7xA0hjgWpKk8XLBqXuAIyTtKGlH4Ii0zMzMKiDLqKo3SDrDX2wpkDS48HV7IqJZ0rkkf/B7AjdExCJJ04CGiKgneTTVF7hVEsCyiJgQEa9J+m+S5AMwLSJeyxC3mZmVUZbEcR3wa0nfAHpI+gjwPySPsDoUEbOAWUVlFxUcH7aJujcAN2SI1czMukiWxPE9ko7xq4H3kPwhvxb4YRfEZWZmVSrLcNwArkx/zMxsC7XJxCHpoIj4Y3p8aHvXRcTscgdmZmbVqaM7jmuAvdPj69u5JoChZYvIzMyq2iYTR0TsXXC8W9eHY2Zm1S7LzPHftFN+e/nCMTOzapdlAuAh7ZQfXIY4zMysm+hwVFU6SQ+SvTimFZ0eCjxf9qjMzKxqlTIct2WdqB5suGZUkKxa++0yx2RmZlWsw8QREZ8DkPRwRPy060MyM7Nq1tE8jiERsTR9+QdJbQ67jYgl5Q7MzMyqU0d3HE8C26bHjSSPp4r3xwiShQvNzGwL0NE8jm0LjrOMwDIzsxrlZGBmZpl01Mcxl3a2eC0UEQeVLSIzM6tqHfVxXJdLFGZm1m101MdxY16BmJlZ99DRo6rPRsQv0uMz2rsu3aHPzMy2AB09qjoF+EV6/Nl2rgm8rauZ2Rajo0dVRxcct7fIoZmZbUGy7DmOpB2AY4D+wArgdxHxelcEZmZm1SnLfhyHAkuBLwH7A/8HWCppXNeEZmZm1SjLHcdVwOSIuKWlQNJJwNXAh8odmJmZVacsM8f7A78uKrsD+ED5wjEzs2qXJXH8HDinqOzstNzMzLYQm0wckuZK+qOkPwL7At+X1CTpUUlNwBXAmFLeSNJ4SU9LapQ0tY3zB0l6XFKzpBOLzq2TtDD9qS/945mZWbllXXJkszZyktSTpC/kcKAJmCepPiIWF1y2DDgdOL+NX/F2RIzenPc2M7PyymvJkQOAxpYNnyTNBCYCrYmjZcMoSevL9J5mZtYFss7jeD9JEtiZgg2dSlhyZADJ/uQtmoADM7z1VpIagGbgkoi4M0NdMzMro5ITh6TjgJuAZ4GRwCJgb+BBOl5ypHjXQChhufYCgyNiRbp17WxJT0bE34rimwxMBhg8eHCGX21mZllkGVX1HeBzETEG+Gf638nA/BLqNgGDCl4PJJl5XpKIWJH+dwkwhzY65CNiekTURURdv379Sv3VZmaWUZbEMTgibi0quxE4tYS684BhknaT1BuYBJQ0OkrSjpL6pMc7Ax+loG/EzMzylSVxvJz2cUCy1MhHgN2Bnh1VjIhm4FzgHuAp4JaIWCRpmqQJAJL2T4f4ngRcK2lRWn0voEHSn4H7Sfo4nDjMzCokS+f4T4GPkcwe/wHJH/H1wPdLqRwRs4BZRWUXFRzPI3mEVVzvYWBUhjjNzKwLlZw4IuJ7Bcc/lzQH2CYinuqKwMzMrDplHY7bExjLu8uqP9IVQZmZWfXKMhx3H+BOYCuSUVIDgXckHR8Rf+6i+MzMrMpk6Ry/gWTZkAERcQDJpL6r8LaxZmZblCyJYzhwZUQEQPrfHwLDuiIwMzOrTlkSxyxgQlHZJ4HflS8cMzOrdpvs45D0C95dGqQnMFPSfJJ1pwYB+wG/6dIIzcysqnTUOd5Y9PovBceLSSb0mZnZFqSjZdX/K69AzMyse8g6j+MQ4LMkI6peAG6KiNldEZiZmVWnkjvHJZ0J3Az8HbgdeBH4laQvdFFsZmZWhbLccVwAHF442U/SzSRrV23WlrJmZtb9ZBmOuxMbL2f+NPC+8oVjZmbVLkvieBC4QtLWAJK2AS4DHu6KwMzMrDplSRxfJFne/A1JLwGvAx8GzuqKwMzMrDqV1MchScB7gcOAD5CujhsRTV0Ym5mZVaGSEkdEhKQngW3TZOGEYWa2hcryqGoByUKHZma2BcsyHHcO8HtJPyNZq6plDSsiwkurm5ltIbIkjo8CzwGfKCoPvCeHmdkWo8PEkQ6//SawGngc+J+IWNPVgZmZWXUqpY/jKpJ9N54CPgVc3qURmZlZVSslcRwFHBERF6THx3ZtSGZmVs1KSRzbRMSLABGxHNi+a0MyM7NqVkrneK90OXW18xovrW5mtuUo5Y7jZZJRU9enP68Wvb6ulDeSNF7S05IaJU1t4/xBkh6X1CzpxKJzp0l6Nv05rZT3MzOzrtHhHUdEDOnsm0jqCVwNHE4y63yepPqIKFxtdxlwOnB+Ud33Ad8C6kiG/s5P667qbFxmZpZdlpnjnXEA0BgRSyJiLTATmFh4QUQsjYgngPVFdY8E7o2I19JkcS8wPo+gzcxsY3kljgEks81bNKVlXV3XzMzKLK/EoTbKoo2yza4rabKkBkkNK1euzBScmZmVLq/E0QQMKng9EFhRzroRMT0i6iKirl+/fpsdqJmZbVpeiWMeMEzSbpJ6A5OA+hLr3gMcIWlHSTsCR6RlZmZWAbkkjohoBs4l+YP/FHBLRCySNE3SBABJ+0tqAk4CrpW0KK37GvDfJMlnHjAtLTMzswrIsjpup0TELGBWUdlFBcfzSB5DtVX3BrwCr5lZVcjrUZWZmdUIJw4zM8vEicPMzDJx4jAzs0ycOMzMLBMnDjMzy8SJw8zMMnHiMDOzTJw4zMwsEycOMzPLxInDzMwyceIwM7NMnDjMzCwTJw4zM8vEicPMzDJx4jAzs0ycOMzMLBMnDjMzy8SJw8zMMnHiMDOzTJw4zMwsEycOMzPLxInDzMwyceIwM7NMnDjMzCyT3BKHpPGSnpbUKGlqG+f7SLo5Pf+opCFp+RBJb0tamP78JK+YzcxsY73yeBNJPYGrgcOBJmCepPqIWFxw2eeBVRGxh6RJwPeAT6fn/hYRo/OI1czMNi2vO44DgMaIWBIRa4GZwMSiayYCN6bHtwHjJCmn+MzMrER5JY4BwPKC101pWZvXREQz8AawU3puN0kLJD0g6eNdHayZmbUvl0dVQFt3DlHiNS8CgyPiVUn7AXdKGhkRb25QWZoMTAYYPHhwGUI2M7O25HXH0QQMKng9EFjR3jWSegHbA69FxJqIeBUgIuYDfwOGF79BREyPiLqIqOvXr18XfAQzM4P8Esc8YJik3ST1BiYB9UXX1AOnpccnArMjIiT1SzvXkTQUGAYsySluMzMrksujqoholnQucA/QE7ghIhZJmgY0REQ9cD3wC0mNwGskyQXgIGCapGZgHfDFiHgtj7jNzGxjefVxEBGzgFlFZRcVHL8DnNRGvV8Dv+7yAM3MrCSeOW5mZpk4cZiZWSZOHGZmlokTh5mZZeLEYWZmmThxmJlZJk4cZmaWiROHmZll4sRhZmaZOHGYmVkmThxmZpaJE4eZmWXixGFmZpk4cZiZWSZOHGZmlokTh5mZZeLEYWZmmThxmJlZJk4cZmaWiROHmZll4sRhZmaZOHGYmVkmThxmZpaJE4eZmWXixGFmZpnkljgkjZf0tKRGSVPbON9H0s3p+UclDSk49/W0/GlJR+YVs5mZbSyXxCGpJ3A1cBQwAjhF0oiiyz4PrIqIPYAfAN9L644AJgEjgfHANenvMzOzCsjrjuMAoDEilkTEWmAmMLHomonAjenxbcA4SUrLZ0bEmoh4DmhMf5+ZmVVAXoljALC84HVTWtbmNRHRDLwB7FRiXTMzy0mvnN5HbZRFideUUhdJk4HJ6cvVkp7OFGF57Qy8sjkV9b0yR1LbNrudwW2dkds6P5Vs611LuSivxNEEDCp4PRBY0c41TZJ6AdsDr5VYl4iYDkwvY8ybTVJDRNRVOo5a53bOj9s6P92hrfN6VDUPGCZpN0m9STq764uuqQdOS49PBGZHRKTlk9JRV7sBw4DHcorbzMyK5HLHERHNks4F7gF6AjdExCJJ04CGiKgHrgd+IamR5E5jUlp3kaRbgMVAM3BORKzLI24zM9uYki/1Vk6SJqePzqwLuZ3z47bOT3doaycOMzPLxEuOmJlZJk4cZmaWiROHmZll4sSRg3TpFMuBpB5u73y4rfMjqWc1tbU7x7uApPeQTGDcKSIqOYO95knaFtgdeE9EzKt0PLXMbZ0PST2AfsBYYHVE/KHovKLCf7idOMpM0knAV0nmnLxOsjzKXOC2dJFGKxNJZ5MsM/MmyZyk7YD7gJ9GxOJKxlZr3Nb5kXQRcDzwKslKGQOAWcBlEdFQydhaOHGUkaQPkkxUPAt4B9gWGALsS5JILomIBRULsIZI6g/8lWT15HdI1jTbEziG5NvaDyPi9spFWDvc1vlJ2/oZ4GPAauCfwD7Amel/r4+ISysXYSKvtaq2FCeRzIS/paVA0jYke5CcAdwq6bCIWFqh+GrJscBjEXF/S4GkR0i+BZ8CXCrpmYj4S6UCrCFu6/wcAjwaEQsLyl6UNBs4Dvi2pLkR8afKhJdw53h5PQv0kbRnS0FE/DMi5kXE2cCjgHcwLI+/ADtKOqqlICLWR8Ty9BvZbDbe88U2j9s6P38G3i/pPyT1aSmMiH9FxK0kyzadXLHoUk4c5XUf8BJwl6TT2xgFMRjYJv+wak9EPEzS3pdL+pqkrWGDEWx7AmsrFV8tSdv6DyRtPdVt3XXSu7argM8A35E0QtJ704EJAPsBL1YswJT7OMpEUs+IWJduazsV+CzQF1gAzAFGAwcCYyLinxULtAYUjiqR9HngApLRPo8ADwF1wA7AxyPirYoFWgNa/l2nx58DzidJFI8CD+K2Lpuif9cnAFNI2rcBeBL4CPAy8MlKt7UTRxlJGhQRy9PjgcAY4JMk+6XfBcyKiCcqGGLNKB6SKOnDJMvxDyH5dvzHiFhSofBqiqRhwAsR8Vb6+GQPkv68YcC9uK27jKRdSFYKfx/J9hTzI+LvlY3KiaMsJO1Bcmt5LDAceAC4HbgnIlak11R87HUtKPoGLKBH8TL7buvykDSCZAju8SR3FQ+S3NX93vM4yit9/Pce4M2Cu46q/XfsxFEGku4kGW57PcleIv8OjCeZBHg1cCmwNiLWVyzIGiHpSmAX4FpgbkubFj4q9H4t5SHptyRbmF4CrAdOAA4GhpJ00l4AvFOtf9y6E0nXkPSB/pLkMeCKiHin4HwfYF1ENFcoxA04cXRS2qfxCjA0IlYVnTsR+Dbw3Yj4ZQXCqylpW79BMsrnwyRJegbwvxGxKL1mOrAsIr5TsUBrQNrWK4C6lsevBecOJunA/Xk1zCno7tK2fpnkbm4YsI6kX7Qe+EtELE//Xb8cEd+sWKAFPKqq87Yh+R/5rOITEXEbcAXwOUk75hxXLfoosBA4MiLeS/KN96PAk5KekHQOyd3eQxWMsVb0AH4H/GfxiYiYA3wZOF5Sv5zjqkVjgaeAMyNiOEmneD9gOnC7pAtJttW+t3Ihbsh3HGUg6RSSP2J3AP8PeCYiXkvP/RswIyJ2rWCINUHS+0hGljxV2BkraWeSiWjTgJXp//mskyQdQfLFZz7J4I55EfF8eu5g4KaIGFi5CGuDpF4kIy5XFC5LJGk7ktn53ybpyxtWmQg35sRRJpI+TfKtoAfJt4dXSZ5Z7gH8KSK+UcHwalLxc19JvwEaI+K8ykZWOyQdRJKUPwD8g6SvY3uSzvL7I2JaBcOrSemjKwoGgfya5PHrRnd/leLEUQbp5Jx/pS+PBA4H+gA7AjcDdxV2dNnmk7QD0CsiXiko60GyfE4jcFhEPFOp+GpF2qYREZE+Zv0YsDewFdAf+BXwYET8axO/xkqQtjVtDZ5JV9qeA3y2moY8O3F0gqSJJBOiXk+LlgG3An8kubWsihEQtaCgrV8jWXF4FcmKob+r9GSoWiNpm/YmqUrq5X/X5VPc1i2z8at9pJoTx2aStB9wJ/Bd4G1ga5IZtXsBzwEXVcNEnVrQTlsPI2nvZcB/tcyXsc5J5yT9iGQrgD+SLG64tuiagRHRVIn4akmJbd06qbiaOHFsJkk/ALaPiDMKynYhWSJgMslyIydExJsVCrFmlNDW2wCfclt3nqQfkazC+nuSkT2rgMeBByLiifSxyq3AJD+m6pwS2lrAbVRhW3tZ9c33N+AoSdtHxBsAEfEyMEvSfJKJPGNJRllZ57it87MTcBnJsi2jSebLfAwYJ+mvJMvo7Fhtf8i6qW7b1k4cm28myWiTmZK+FRGPtZyIiJck7YpXwi0Xt3UO0lFqtwOr0gmViyTdTtIpvj+wK8nAjwmVi7I2dPe29qOqzdDSQShpN5Ln7geTzB6fRTJBbRzJaqGeT9BJbuv8SeoTEWvaWEjyIGBORHjicJl017Z24thMkrZreaYuaW/gE8CngJ2BW0gWgquK/YG7O7d1fiRtGxH/KCprSd5nAPtFxDkVCq+mdOe2duLISNJeJKvgTiKZBHUXydIMcz1Xo7zc1vkpauvtgLuB3wIPtQwXTbcKWB0Rr7f7i6xDtdDWThwZSZpDssvfzcC2JNs4jknLpkXEHV6htTzc1vlpo61PJNlt7u8kQ8t/W7noaksttLUTRwbpENAlEdG3qHw7kslpk4AzIuLBSsRXS9zW+XFb56dW2roqO16q2Fbmwpq0AAABtklEQVTAPEnHFhZGxJsRcRHJmOsz02UCrHPc1vlxW+enJtraiSOb5SRLdn9L0n9IGinpvQXnFwEjqnHcdTfkts6P2zo/NdHWnseRQbrg27dIVgg9CPgQ8KKkt0lWDD0auKmCIdYMt3V+3Nb5qZW2dh9HiSQNJ1neYnuSO7XdSVbAbSKZVzASuAa41Z21neO2zo/bOj+11NZOHCWS9BTJYmRvkqwpsyMwCFgDTK/2zqzuxG2dH7d1fmqprZ04SiDpSODqiNgjfd0LGEAyhO4Ykv/xT/cKrZ3nts6P2zo/tdbW7hwvzTbAS5IGAUREc0Q8HxG3AxeS7A9xZCUDrCFu6/y4rfNTU23txFGaO0gm5/xI0tDCE+k3hL+RrM5qnee2zo/bOj811dZOHCVIFx/7vySj0BZKmivpy5JGSTqXZAXL/61okDXCbZ0ft3V+aq2t3ceRkaR9gYnACcAHgdkki+zdUNHAapDbOj9u6/zUQls7cXRCOnGnd8vmQtZ13Nb5cVvnp7u2tROHmZll4j4OMzPLxInDzMwyceIwM7NMnDjMzCwTJw4zM8vEicPMzDJx4jAzs0z+P/aZ2vmjXlE0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2271066ee10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "q = QuantumRegister(5, 'q')\n",
    "c = ClassicalRegister(5, 'c')\n",
    "\n",
    "shor = QuantumCircuit(q, c)\n",
    "circuit_aperiod15(shor,q,c,7)\n",
    "\n",
    "backend = Aer.get_backend('qasm_simulator')\n",
    "sim_job = execute([shor], backend)\n",
    "sim_result = sim_job.result()\n",
    "sim_data = sim_result.get_counts(shor) \n",
    "plot_histogram(sim_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see here that the period, $r = 4$, and thus calculate the factors $p = \\text{gcd}(a^{r/2}+1,15) = 3$ and $q = \\text{gcd}(a^{r/2}-1,15) = 5$. Why don't you try seeing what you get for $a = 2, 8, 11, 13$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Free flow\n",
    "You could now try to implement this algorithm for large integers and even conceive a code structure that would handle any integer (python is good because you don't have to make integer memory size considerations), and you are encouraged to do so.\n",
    "\n",
    "The bad news is your code would not be able to run in a quantum simulator (which can simulate quantum properties at an exponential cost) nor in a real quantum device on IBMQX (as the number of qubits would be insufficient).\n",
    "\n",
    "However, you may write this code and leave it as a legacy for your grandsons and maybe, just maybe, a time will come when they execute it and it works well enough to break nowadays encryption mechanisms like RSA."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
